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Abstract 

The wavelet transform modulus maxima (WTMM) used in the singularity analysis 
of one fractal function is extended to study the fractal correlation of two multifractal 
functions. The technique is developed in the framework of joint partition function 
analysis (JPFA) proposed by Meneveau et al. [1] and is shown to be equally effective. 
In addition, we show that another leading approach developed for the same purpose, 
namely, relative multifractal analysis, can be considered as a special case of JPFA at 
a particular parameter setting. 

1 Introdcution 

Fluctuations in many natural and artificial phenomena are found to exhibit fractal charac- 
teristics. In applications, this has been characterized by the so-called singularity spectrum of 
some numerical or experimental data. To understand the dynamics underlying the fractal, 
it is not uncommon that multiple data capturing different aspects of the phenomenon of 
interest are used in the analysis. For example, velocity and temperature flucutations are 
used to analyze the momentum and energy aspects of the multifractal hydrodynamic tur- 
bulence [1], blood pressure and heart rate fluctuations to analyze the cardiovascular aspect 
of the 1 / /-like power spectrum of the heart rate variability in humans [2] , packet size and 
arrival time to analyze the congestion and connectivity aspects of the multifractal network 
traffic [3], and so on. For fractal analysis on natural objects, multiple data cross examination 
may provide the chance to examine the potential link of fractal fluctuation in the data. In 
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particular, one would suspect some degree of fractal correlation in the data if the fractal 
generating mechanisms associated with the data source are coupled together. 

Essential to the notion of fractal correlation is the distinguishability of singularity spectra. 
There are fundamental and practical issues related to this subject. For example, consider 
a standard iVp-adic multinomial process on an interval. It is a multiplicative cascade con- 
structed by repeatedly dividing the interval into equal N p segments and assigning (probabil- 
ity) weights Pi,i = 1, ■ • • , N p , from one generation to the next. Continuing this procedure ad 
infinitum leads to a limiting process with no density (almost surely) and intermittent spik- 
ing pattern. Its singularity spectrum /^(a) may be estimated by the Legendre transform of 
r w (q) = -\og(Epf)/log(N p y. 

f w (a) = qa- r n {q) 

where a(q) = dr^/dq. Now, consider a different iV m -adic cascade (N m ^ N p ) generated by 
weights mi, j = 1, • • • , N m , and its singularity spectrum / M (a). If f n (a) = / M (a), one must 
have 

However, no {pi} and {mj can be found to satisfy this equation for all q. Thus, singular- 
ity spectra can in theory be distinguished, at least for the important class of multinomial 
processes. For more in-depth treatments and examples, see the excellent book by Pesin [4]. 

In practice, a different issue can arise. That is, two singularity spectra may be close 
to each other within the limit of finite precision. Consider again the cascades from above. 
Let 7v(-oo) = r M (-oo), r^+oo) = r M (+oo) (so max({pj) = max^m,}) 10 ^/ 10 ^™), 
min({pj}) = min({mi}) log ^ JVp ^ log ^ Arm - ) ). Then, f w (a), f^a) will agree at four important q 
values: — oo, +oo, 0, 1. With the rest of pi and m« chosen properly, they can be made 
almost indistinguishable (Fig. 1). This problem was addressed by Levy-Lehel and Vojak who 
developed a much sharper mutual multifractal analysis to relate the singularity spectra to the 
generation of multinomial processes [5]. Riedi and Scheming arrived at the similar relative 
multifractal analyses with further details on the numerical implementation [6]; see also [7, 8]. 
For experimental data, Meneveau et al. introduced a joint partition function analysis (JPFA) 
based on the ID version introduced by Hentschel and Procaccia [9] and Halsey et al. [10]. 
With essentially the same procedure of estimating the singularity spectrum, these authors 
characterized fractal correlation in the small scale kinetic energy transfer, heat concentration 
and vorticity of the turbulent flow [1]. They implied that multiple cascades of more than 
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one variables are responsible for the fractal fluctuation in fluid turbulence; see also [11]. The 
JPFA has also been applied in diversed areas such as precision agriculture [12], soil property 
[13], and re-emerged in the general discussion of discrete scale invariance of the multinomial 
process [14]. 

In these past studies, JPFA as well as most singularity analyses were conducted on 
the assumption of one-dimensional multifractal measures and solved using the classical box 
counting procedure. Such applications can be limited in scope when, typically in the experi- 
mental study, the data is a fractal function that is not additive. For functions created by the 
integral of multinomial measure under C°° perturbation [15], Bacry et al. proved that the 
singularity analysis developed for fractal measures is equally applicable by using the so-called 
wavelet transform modulus maxima (WTMM) method. Later, the validity of WTMM was 
examined by Jaffard for any function [16]. Regarding WTMM, it was proved that (i) it can 
yield the upper bound estimate of the singularity spectrum for any function and (ii) it is 
exact for the so-called self-similar fractal function, as long as the so-called maxima lines are 
not too close to each others. 

The purpose of this study is to introduce a joint WTMM method to carry out the JPFA 
of fractal correlation. The WTMM-based approach can provide an accessible and stable 
tool for estimating fractal fluctuation in multiple experimental data. We also point out that 
JPFA is a more general formulation. In particular, the existing relative multifractal analysis 
developed for the similar purpose relates to JPFA at a particular parameter setting. While 
the term correlation is normally linked to the second order statistics, fractal correlation 
as estimated from the singularity spectra is a property of moment of all orders. Indeed, 
the primary object of the analysis is the Hausdorff dimension / (011,012) of the support of 
observing Holder exponents Oi and 02- It will be shown that the set /(oi,02) describes 
a two-dimensional surface whose level sets characterize the coupling of fractal generating 
mechanisms as well as the relative multifractal spectrum studied in the past. 

Our results are organized into four sections. In the next section, the background of 
WTMM is first summarized. The extension to the WTMM-based JPFA and its connection to 
the relative multifractal spectrum are then given. The test of the method using multinomial 
cascades are presented in section 3. Concluding remarks are given in section 4. 
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2 WTMM-Based Fractal Correlation Analysis 
2.1 WTMM Singularity Analysis 

Singularity analysis is built on the notion of Holder continuity of functions. Recall that a 
function x(t) is Holder continuous of exponent a' if there are a', Sq, C G R + , such that, for 
S < 6 , 

\x(t + S)-x(t )\ <C\S\ a '. 

In the neighborhood of x(to), there exists a supremum a(to) that (1) is valid for all a' < a(t ). 
The exponent a(t ) is the Holder exponent of x(t) at to- Formally [15, 16], one can find an 
nth order polynomial P n {t) and a(t ) G [n,n+ 1) such that 

\x{t + S)-P n {t )\<C\S\ a ^. (1) 

It is evident that the Holder exponent characterizes the differentiability of the function 
and, thus, how the function can fluctuate. For example, a = +oo for C°° functions, a G 
(n,n + 1) for functions that are only n times differentiate and a < 1 for functions that are 
non-differentiable. The a < 1 case draws the most attention since it means the function 
can fluctuate in large amplitude over short time intervals and gives rise to the so-called 
intermittent pattern witnessed in many physical systems. In this case, the Holder exponent 
a(t) is also known as the singularity exponent. 

The natural tool to analyze the singularity property is by the wavelet transform: 

Tj,[x](t,a) = - r V (—) x{t')dt' (2) 

CL J —oo V CL J 

where T^[x](t, a) is the wavelet coefficient and ip(t) is the analyzing wavelet. Muzy et al. 
showed that the exponent a(t) can be estimated effectively using the supremum of |T^,[x]| 
along the so-called maxima line formed by the local wavelet modulus maxima [18]. Denote 
the set of maxima lines at scale a by C(a) = {li, I2, ■ ■ ■ , Zjv(o)}- Bacry et al. proved that [15] 

Z{a-q)= £ C?~a^ (3) 
keC(a) 

where Cj = swpu a \ el . \ T^[x] (t, a) \ is the supremum of the modulus maxima along the maxima 
line For functions created by the integral of multinomial measures under C°° perturba- 
tions, it was shown that the Legendre transform of r(q) yields the Hausdorff dimension of 
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the support {t, a(t) = a}, /(a); 

r(q) = mm(qa - /(a)); (4) 

see also [16]. In the literature, Z(a; q) is sometimes referred to as the partition function due 
to its analogy to the energy partition function in statistical mechanics. A monofractal refers 
to the case where {a} is a singleton. The singularity spectrum is called multifractal when 
{a} spans an interval. 

2.2 WTMM-based JPFA of Fractal Correlation 

In this section, the WTMM approach is generalized and applied in the JPFA of fractal 
correlation between multifractal singularity spectra. We present the application of JPFA of 
two data sets. The extention to more data sets is conceptually similar. 

Consider xi(t),x 2 (t) and their respective sets of singularity exponent {«i}, {a 2 }. Let the 
maxima lines of | at scale a be denoted as Ck{a), k = 1, 2. A natural extension of the 

existing WTMM analysis is to consider a joint partition function of the form: 

Z(a;q 1 ,q 2 ) = Y,CZ ij fZ u) (5) 
j 

where C^ r) C 2 s are the modulus maxima along the maxima lines G £i, l 2tS £ £2- 

To realize (5), the maxima lines in £fc,fc = 1,2 must be paired up properly (so the 
index j can run). As in most correlation analyses, the goal is to characterize the property 
related to observing both singularity exponents a\ and In terms of the WTMM analysis, 
such information is contained in the modulus of the neighboring maxima lines. If the time 
coordinate of /fcj(a) is denoted by t k j(a), this means the coefficients Ci >r , C 2tS paired up in 
(5) can be determined by 

\ti, r ~ t 2 . s \ = mm(\t iy - t 2 ,s\) = min(|ti ir - t 2 , s '\) (6) 
r s' 

Once (5) and (6) are established, similar procedure developed by Bacry et al. can be 
extended to characterize the geometry associated with the observation of a\ and a 2 . In 
particular, based on Ck : \ ~ a ak ^ x \ A = r, s [15,17,18], (5) can be given by 

Z{a; qi ,q 2 ) ~Y. a9iai{rij))+92a2(s{j)) = J J da 1 da 2 V{a 1 ,a 2 )a qiai+g2a2 a - f{ai ' a2) (7) 
j 
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where V(aci, a 2 ) and f(a\, ot 2 ) are the probability density function and Hausdorff dimension 
of the support of (ai, a 2 ), respectively Applying the standard argument of steepest descent 
in small a, one has 

Z(a; qi ,q 2 )~a T ^ (8) 

where 

r (<?i, Q2) = min(giai + q 2 a 2 - f(a u a 2 )). (9) 
Hence, r(qi,q 2 ) and f(ai,a 2 ) are Legendre transform pair: 

at = dr(q 1 ,q 2 )/dq 1 ,a 2 = dr(q 1 ,q 2 )/dq 2 , 
/(ai, a 2 ) = ax(qt, q 2 )q x + a 2 (<?i, q 2 )q 2 - r(<?i, q 2 ). (10) 

where 

<?i = df/dai, q 2 = df/da 2 . (11) 
Finally, from (9) ~ (H), the correlation coefficient between a x ,a 2 can be estimated using 
7"(gi,? 2 ): 

cov(ai,a 2 ) 



P 



2 r d 2 r 



(12) 

91=92=0 



where cov denotes the covariance and a\ denotes the standard deviation of A. This expression 
will be used in the next section to compare with the numerical result. 

In practice, the Legendre transform (10) relies on using r estimated from (8). However, 
there are known factors, such as lacunarity [19, 20], that introduce oscillatory, scale depen- 
dent, prefactor. This results in the poor estimate of r(qi, q 2 ). A remedy to this problem can 
be motivated by an alternative approach equivalent to the canonical ensemble in statistical 
mechanics [20]. Let 

fin fii'i 
Z{a; q 1 ,q 2 ) 

Then, it can be shown (Appendix): 

Ai(a; ft, q 2 ) = E Ki, «! Qi, ft) Md,^)) ~ a a ^^ (14) 
3 

A 2 (a; qi, q 2 ) = E "0\ "I 9i, ft) log(C 2 , s(j) ) ~ a Q2 ^ 2 ) (15) 
i 

F(a; gi, q 2 ) = J2 a ; ft> ft) log(^(j, a; ft, 92)) ~ a /(ai,a2) . (16) 



6 



2.3 JPFA and Relative Multifractal Analysis 

Relative multifractal analysis and similar ideas were developed to characterize fractal corre- 
lation between fractal measures. The main idea is to replace the use of Lebesgue measure 
in the traditional fractal analysis [5,6,7,8]. Specifically, consider the partition functions of 
multifractal measures tt and /x 

J2 n(A) q ~ \A\ r ^ q \ ^ A ') q ~ I^T m(9) ( 17 ) 

A£H A'eW 

where H, H' denote generic partitions of the support and | ■ | denote the Lebesgue measure 
of the set. To examine the extent to which the singularity of tt correlate with /x, the sets 
which scale as a power law will now be characterized by using tt. For example, the partition 
function of tt is now written as 

Z*(A)WA)-*® ~ 0(\A\) (18) 

where the "big O" describes the order relationship 0(|A|) — > const, as \A\ — > 0. Define 
T n /^(q) = sup{t(g)} for which (18) holds. The relative multifractal spectrum is obtained via 
the Legendre transform of Tw^. It characterizes the support of the singular behaviour of the 
form it ~ y^T/M) where a^/^q) = dr 7T /^(q)/dq. The relative multifractal analysis can draw 
a much sharper distinction between it and /x. For example, T n /^{q) is nonlinear when tt ^ /x 
and T n /fj,(q) = q — 1 when tt = /x; see [6] for more details. 

Comparison of (18) with (5) and (8) suggests T n /^ can be obtained as the level set of 
r (9i) Q2) = where 

<?2 = -tV/m^O- ( 19 ) 

To assure finite generalized dimension, tv/^I) = and, thus, the level set r(qi, q?) = must 
pass through (gi,^) = (1>0). This property can also be directly seen for the multinomial 
processes (next section). 

Similarly, by switching the role of qi,q%, the singular behaviour of /x "viewed" by tt can be 
described. With the same argument, this is characterized by 7^/^(92) defined from the same 
level set r(qi, 52) = where q\ = — r M / 7r (g 2 )- It may be useful to point out that 7V/ M , r^/^ on 
the qi x q 2 plane are nothing but mirror images of the contour of r(gi, q%) — about q 2 = 
and qi = axes, respectively. In general, r n /^(q) 7^ T /J / 7r (g), although they are derived from 
the same level set. 
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3 Numerical Experiments 



To test if the WTMM-based JPFA can reliably characterize fractal correlation, numerical 
experiments are conducted on the coupled random binomial cascades studied by Meneveau 
et al. [1]. 

The first cascade, denoted as ir, is generated by weights po,p± (referred to as 7r-cascade). 
Let I ri ,—,rj denote an interval segment generated in the Jth iteration where G {0, 1} and 
t{1~ % is the based-2 coarse-grained representation of any x G I rit ... trj . By the multiplicative 
rule, ir(I rii ... irj ) = Y[j=iP rj - The second cascade, denoted as /i, is generated by weights m , mi 
(referred to as the /i-cascade). With the same addressing scheme, one has fJ>{I sl ,—,sj) — 
ll/=i m Sj where Sj G {0, 1}. 

In the numerical experiment, a parameter g and a uniform random variable 7 in [0,1] are 
used to control the degree of coupling or correlation between the cascades. Let 1^, Ir be 
the new segments created from their parent segment of the previous generation. If 7 < g, 
the weights assigned to Il,Ir of the //-cascade will depend on exactly how the weights of 
the 7r-cascade are assigned. The rule for this dependence is that po and m (pi and mi) are 
always assigned at the same time. For example, if po is assigned to II (Ir) of the 7r-cascade, 
mo will be assigned to II (Ir) of the //-cascade and similarly for p\ and mi. If 7 > g, the 
weight assignment for the cascades will be completely independent from each other. This 
way, the fractal generating mechanisms of the cascades are completely dependent of each 
other when g = 1 + and independent of each other when g = 



The coarse-grained joint partition function for the coupled cascades can be defined based 
on (5): 



From the combination of 7 completely dependent and (1 — 7) independent proportions, Zj 
is derived explicitly as 



Zj{a;q u q 2 ) = XMW-,o) 9l M-fn,---,o) 



Zj(a- qi ,q 2 )~(2Y) 



.1 



(20) 



where 



As J 



Y = 7 



00, Z, 




(21) 



r(qi,q 2 ) 



log 2 (2F) 



(22) 



S 



By (10), the analytical a±, a 2 and f(cei, a 2 ) can be found [21]. These results will be compared 
to the numerical ones below. 

In the numerical experiments, po = 0.2, p\ = 0.8 and mo = 0.4, mi = 0.6 are used to 
generate the tt- and /i-cascades for g = 1,0.8,0.3,0. For each g value, 30 pairs of n, fi 
cascades, each of 16,384 points are generated. The first derivative of the Gaussian wavelet 
has been used as the analyzing wavelet in this work. Using higher order derivative of the 
Gaussian wavelet does not create qualitatively different result. In practice, the modulus 
maxima and the maxima lines of the individual cascade are first obtained. The modulus 
maxima from the nearest maxima lines are then paired up according to (6) and used in 
(5) to define the joint partition function. The numerical ai,a 2 and /(tti,a 2 ) ar e finally 
estimated following (13) ~ (16) and r(gi,g 2 ) is determined following (9). 

Typical maxima lines of the coupled cascades are shown in Fig. 2. It is observed that the 
maxima lines are "aligned" when the fractal generation is completely dependent at g — 1 
and begin to "mis-align" for g < 1. The power law scaling of Z(a;gi,g 2 ) are found in all 
cases (Fig. 3). In Fig. 4, contours of the level set of /(ai,a 2 ) are shown on the ot\ x a 2 
plane. It is evident that the geometry of the contour lines vary systematically with the 
g value. When the fractal generating mechanisms are completely dependent of each other 
(g — 1), /(cm, a 2 ) describes a one-dimensional curve supported by the functional relationship 
ai(a 2 ). This is expected as any spiking pattern in one cascade automatically implies the 
same for the other. As a result, the maxima lines will converge at the same location in the 
time-scale plane. This establishes the one-to-one relationship of observing the exponents a\ 
and a 2 . In general (g < 1), /(ai,a 2 ) describes a two-dimensional surface, which gives rise 
to oval-shape contours (Figs. 4b, 4c, 4d). This means that the observation of a\ can take 
place simultaneously for a range of a 2 . As a result, the contour "opens up" and becomes the 
largest when the fractal generations are completely independent from each other {g = 0). 
For g = 0, there is a perfect alignment of the axes of the contour and the a\ = and 
a 2 = axes (Fig. 4d). Superimposed on these figures are the analytical /(ai,a 2 ) derived 
by the Legendre transform of (22). It is seen that the WTMM-based JPFA agrees well with 
the theory. To further the check of the WTMM-based approach, the correlation of ai,a 2 
is estimated using the numerical T(gi,g 2 ). This is to compare with the analytical p = g 
obtained by substituing (22) into (13). The result is summarized in Table 1. Again, good 
agreement is found. 
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9 


0.0 


0.3 


0.8 


1.0 


P 


0.016 


0.280 


0.775 


0.838 



TABLE 1 Numerical p value for g — 0, 0.3, 0.8, 1.0 (Note, in theory, p = g). 

To test the robustness of the algorithm, different analyzing wavelets are also used to 
study the coupled cascades. While deviations are expected to result from the choice of the 
analyzing wavelets, no qualitatively different result is found. Fig. 5 demonstrates the coupled 
cascades of g = 0.8. It is evident that the contour lines of f(a±, a 2 ) estimated from different 
analyizing wavelets are all failed onto the theoretical contour lines. 

Finally, from the numerical data, 1^/^(91) is estimated from the level set r(qi, q 2 ) = 0. We 
then conduct the relative multifractal analysis using the deterministic algorithm proposed by 
Riedi and Scheming [6]. In Fig. 6, the r^/ M (gi) estimated by these two different approaches 
are shown to match well. The agreement confirms that the relative multifractal spectrum is 
contained in the level sets of r(q l7 q 2 ). Similar match is also found for t^/^ (not shown). 

4 Concluding Remarks 

In applications, the ability to characterize fractal correlation in the data fluctuation could 
provide insights into the underlying complex dynamics. 

In this work, a WTMM-based technique is introduced for the first time to estimate the 
fractal correlation in the framework of joint partition function analysis proposed by Meneveau 
et al. [1]. As WTMM has been proven an effective tool to extract the singularity spectrum of 
certain important class of fractal functions [15, 16], it is shown that the extension developed 
in this work also capture accurately the fractal correlation of data fluctuation. We also 
found another leading idea developed for the fractal correlation analysis, relative multifractal 
spectrum, can be considered as a special case of JPFA at a particular parameter setting. 

Acknowledgment This research is supported by Natural Science and Engineering Research 

Council of Canada. 
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First, consider (14). By (8) and (10), one has 

dZ(a; qi , q 2 )/d qi ~ a r(<?1 ' 92) log(a)<9r(gi, te)/d qi = a T ^^ log(a)ai. (Al) 
Note the prefactor log(a)ai in (Al) that varies logarthmically with a. From (Al), one has 

E-MI^I) = ^. (A2) 

Equation (14) follows after substituting (Al) into (A2). Note also the prefactor in (8) is 
canceled out in (A2). The derivation for (15) is the same. 
Based on (A2), one has 

E ,lo gM = ?1 ^l + ?2 ^*-log(Z). (A3) 
Again, using (Al), (16) results after substituting (10) into the above. 
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Figure Captions 



Fig. 1. Multifractal analyses of 3-adic 7r-cascade of weights 0.2, 0.28, 0.52 (symbol "o") and 
5-adic //-cascade of weights 0.09463, 0.1200, 0.1800, 0.2217, 0.3837 (symbol "+"). (a) r n {q) 
and 7^(9); (b) f n (a) and f^a). 

Fig. 2. Typical maxima lines in the time-scale plane from one of the 30 sets of it- ("o") 
and //-cascades ("+") with coupling parameter g = 1.0,0.8,0.3,0.0 (top to bottom). Notice 
the perfect alignment of maxima lines for the completely dependent cascades (g = 1). First 
derivative of the Gaussian wavelet is used in the numerical calculation. 

Fig. 3 log(Ai), log(Ai) and log(F) vs. log(a) plots of a typical case of the coupled cascades 
with g = 0.8. The straight lines describes the power laws at (91,92) = (3,-2), (4,0), (0,0), 
(-1,3) (top to bottom). Regression lines are shown as solid lines. The slope of the regression 
lines are estimated as cci, «2 and /(ai,^) based on (14)-(16). 

Fig. 4 Averaged contour lines of numerical /(ai,a 2 ) — C for C = 0.6 (outer contour), 0.75, 
0.9 (inner contour). Theoretical contour lines are shown as solid lines. The averaging over 
30 pairs of tt and /1 cascades is shown as "o" (C = 0.6), "•" (C = 0.75) and "A" (C = 0.9). 
Error bars of one standard deviation from selected data points are shown. They are obtained 
from the ensemble of contour line points in a uniform grid on the 9 x p. The corresponding 
g values are (a) 1, (b) 0.8, (c) 0.3, (d) 0. 

Fig. 5 Averaged contour lines of numerical f (011,0.2) = C for g = 0.8 and C = 0.6 (outer 
contour), 0.75, 0.9 (inner contour). The analyzing wavelets are ... ("o"), ... ("+"), ... ("A"). 
Theoretical contour lines are shown as solid lines. Error bars of one standard deviation from 
selected data points are shown. Compare with Fig. 4b. 

Fig. 6 7^/^(9) estimated by WTMM-based JPFA (V) method (from the contour line of 
r(q,p) = 0) and the deterministic algorithm porposed (solid line); see text. The solid lines 
shown are based on the ensemble average with ± 3 standard deviation boundaries plotted 
as long-dashed lines. The g values are (a) 1, (b) 0.8, (c) 0.3, (d) 0. 
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